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We explore and describe the roles of inter-molecular vibrations employing a Brow¬ 
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bath interactions, which we use to analyze two-dimensional (2D) THz-Raman spec¬ 
tra obtained by means of molecular dynamics (MD) simulations. In addition to 
linear absorption (ID IR), we calculated 2D Raman-THz-THz, THz-Raman-THz, 
and THz-THz-Raman signals for liquid formamide, water, and methanol using an 
equilibrium non-equilibrium hybrid MD simulation. The calculated ID IR and 2D 
THz-Raman signals are compared with results obtained from the LL-I-SL BO model 
applied through use of hierarchal Fokker-Planck equations with non-perturbative and 
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the signals obtained from the MD simulations are reproduced with the LL-I-SL BO 
model, indicating that this model captures the essential features of the inter-molecular 
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izability, and dephasing time. The origins of the echo peaks of the librational motion 
and the elongated peaks parallel to the probe direction are elucidated using optical 
Liouville paths. 
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I. INTRODUCTION 


Molecular vibrations in condensed phases play an essential role in various dynamic pro¬ 
cesses, including inter- and intra-molecular couplings, and solvent dynamics, all of which 
entail energy exchange as well as thermal excitations and relaxations.^ Multidimensional 
vibrational spectroscopy techniques make it possible to experimentally distinguish such pro¬ 
cesses due to the sensitivity of the nonlinear response functions utilized in these techniques 
to complex dynamics.® For intra-molecular vibrations, the roles of relaxation and dephas¬ 
ing are well understood both theoretically and experimentally due to the advent of infrared 
(IR) laser technologies. Methods of analysis with theoretical models that utilize molecular 
dynamic simulations have also been developed to elucidate multidimensional IR signals.® 
Because the primary inter-molecular modes, which are the objects of study in 2D IR spec¬ 
troscopy, can be separated from the other modes, as in the case of the OH stretching mode 
in liquid water, stochastic models whose parameters are obtained from classical molecular 
dynamics (MD) simulations have been useful for analysis of the inter-molecular vibrational 
modes. For inter-molecular vibrational modes, two-dimensional Raman spectroscopj^ was 
for a long time the only two-dimensional spectroscopy that could be used for experimental 
study. However, due to technical difficulties, such investigations have been carried out only 
for CS 2 ,^^® Benzene,!^ and formamidd^ liquids. Theoretical investigations have also been 
limited, due to the availability of experimental data and limitations on computational power 
for simulations. 

Two-dimensional THz-Raman spectroscopy, which has been studied both theoreticalljh^^ 
and experimentally,^ has created a new possibility for investigating the details of inter- 
molecular vibrations. In 2D Raman spectroscopy, the observable is dehned in terms of 
the three-body response function for the polarizability of the system, IT -^RRR(^2,ti) — 

— ([[11(^2 + H), n(ti)], II(0)])/h^, where (...) represents the thermal average and A(t) = 

is the Heisenberg operator for an arbitrary operator AW In the case of 
2D THz-Raman spectroscopy, the response function consists of one polarizability, II, 
and two dipole operators, /i, and there are three different measurements, which de- 
pend upon the sequence of the Raman and THz pulses as = —{[[fi{pt 2 + 

H),/i(H)],n(0)])/h2, = -{[[fi{t 2 + ti),n{ti)],fi{0)])/h^, and Rttr(^ 2 ,H) = 

— ([[11(^2 + H), A(H)], A(0)])/h^. While inter-molecular vibrational modes are usually both 
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Raman and IR active, the types of information that we can obtain from the 2D Raman 
signal and each of three THz-Raman signals are different, due to the role of the nonlinear 
polarizability. Because each of the above-mentioned response functions is dehned in terms 
of the three-body correlation function, the signal will vanish if the system is harmonic and if 
the total dipole moment and polarizability are linear functions of the collective coordinate, 
q, representing the inter-molecular vibration, because there is an odd number of Gaussian 
integrals involved in the response function; Tr{g(t2 + ti)q{ti)q{0) exp{—/3Hs)}. The dipole 
moment is approximated reasonably well as a linear function of q as fi{q) = Hiq, because the 
total dipole moment is a linear function of the distance between the charges in the system, 
and the nonlinear dipole-induced dipole interactions are weak. However, the contribution 
of the non-linear polarization is not negligible, because the polarizability originates in the 
electronic states of molecules, which depend on the complex conhgurations of the atoms 
and molecules. For this reason the polarizability is expressed in a Taylor expansion form as 
n(g) = Y\.iq + I\. 2 q ^Because n has this non-linear form, the three response functions give 
above, representing the observables in 2D THz-Raman spectroscopy experiments, provide 
information about three different physical processes.^ Contrastingly, because 2D Raman 
spectroscopy experiments measure just a single observables, they do not provide such a 
detailed picture of the physical system. The richness of the information obtained through 
2D THz-Raman spectroscopy allows for a detailed analysis of inter-molecular vibrational 
modes. Note that the optical setup for the TTR measurement differs signihcantly from those 
for the RTT and TRT measurements. This is because the RTT and TRT responses are 
detected as the emission of THz signals, while the TTR response is detected as an induced 
Raman signal. 

Although we can obtain relatively reliable 2D THz-Raman signals using the full MD 
simulation techniques developed for 2D Raman spectroscopy,^^^^ analysis of the spectra is 
not straightforward, due to the complexity of the 2D prohles of the signals, which arises 
from the complexity of the inter-molecular vibrational modes. As demonstrated by 2D 
Raman and 2D IR spectroscopy studies, a model-based analysis is useful for treating this 
problem, because the 2D prohle of the signal is so sensitive to the underlying dynamics 
that the complex 2D prohle cannot be reproduced without capturing the essential features 
of the vibrational modes.^^^ While stochastic models, which can be regarded as Brownian 
models with non-linear system-bath interactions,^^^® are recognized as versatile models for 
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analyzing intra-molecnlar modes observed in 2D IR spectroscopy experiments, it is not clear 
if snch models are nsefnl in the 2D THz-Raman case. This is becanse in contrast to the 
intra-molecnlar modes, which are clearly definable with the normal mode pictnre, the inter- 
molecular modes are not localized and change in time due to changes in the configuration 
of the system molecules. 

In this paper, we explore the possibility of characterizing inter-molecular modes using 
a Brownian model with linear-linear (LL) and square-liner (SL) interactions utilizing 2D 
THz-Raman signals obtained from MD simulations. In order to treat a non-perturbative, 
non-Markovian, and nonlinear system-bath interaction, which is necessary to describe the 
effects of homogeneous and inhomogeneous broadening in a unihed manner, we employ the 
hierarchal equations of motion approach.®® The properties of inter-molecular motion are 
investigated using the htted model. 

This paper is organized as follows. In Sec. |TT| we explain the methodology for calculating 
2D THz-Raman signals from full MD simulations. In Sec. H we present the LL-I-SL BO 
model and the hierarchal equations of motion formalism. We then show how this formalism 
can be used to calculate 2D signals. The MD and fitted results obtained from the LL-I-SL BO 
model are presented and analyzed in Sec. W Section |V] is devoted to concluding remarks. 


II. FULL MD SIMULATION 

While, to this time, the experimentally obtained 2D THz-Raman signals are limited to the 
case of liquid water, in this paper we analyze 2D signals obtained from full MD simulations 
for formamide, water, and methanol. We chose these liquids from among many substances 
that have been investigated in full MD studies of the 2D Raman and 2D THz-Raman spec¬ 
troscopy as characteristic examples of 2D THz-Raman signals. The MD simulation results 
used in the present study of these molecules for 2D Raman-THz-THz (RTT) and THz- 
Raman-THz (TRT) signals were originally presented in a previous study.® Nevertheless, 
here we repeated the full MD simulations in order to also obtain THz-THz-Raman (TTR) 
and infrared absorption signals, in addition to the RTT and TRT signals. Moreover, we 
employed the Ewald sum for the evaluation of the dipole and polarizability, in addition to 
the force fields. This contrasts with the situation in previous studies, in which only force 
helds were computed with the Ewald sum. The change in the resulting signals due to the use 
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of the Ewald sum in the computation of the dipole and polarizability, however, are small. 


A. Models and simulation details 

Based on the MD simulations, we calculated the linear absorption (IDIR) spectrum and 
2D THz-Raman signals of liquid formamide, water, and methanol. Each system consisted of 
108 molecules in a cubic box with periodic boundary conditions. The interactions between 
the molecules were modeled by a modihed T potentialj ^^ * ^° l the TIP4P/2005 potential,^ and 
the B3 potentiaP^ for formamide, water, and methanol, respectively. 

The interaction potentials were cut off smoothly at a distance equal to a half the length 
of the system using a switching function, and the long-range Coulomb interactions were 
calculated with the Ewald sum. The intra-molecular geometries were kept rigid throughout 
the simulations, using a constraint provided by the RATTLE algorithm. The equations 
of motion were integrated using the velocity-Verlet algorithm with time steps of 5.0 fs for 
formamide and 2.5 fs for water and methanol. The system volume and total energy were 
hxed after the completion of the isothermal simulations carried out for equilibration. The 
conditions of the simulation were set such that the average densities were 1.120 g/cm^ for 
formamide, 0.997 g/cm^ for water, and 0.786 g/cm^ for methanol. The temperature was set 
to 300 K. The permanent molecular polarizability of each liquid was utilized with the atomic 
polarizability for formamide and methanol,®! and the Huiszoon polarizability for water.®! 


B. Molecular polarizability and dipole 

While the MD simulations were carried out using the permanent polarizability, we calcu¬ 
lated the 2D THz-Raman signals using a full-order dipole-induced-dipole (DID) polarizabil¬ 
ity model. We did this because the 2D prohles are extremely sensitive to the accuracy of the 
calculated optical observables. In the DID polarizability model, the expression determining 
the polarizability of a molecule includes contributions from other molecules, and interactions 
between molecules are dehned with respect to the centers of individual molecules.®®! The 
total polarizability of the system in a MD simulation is given by n(t) = Ilj, where 11 j 
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is the polarizability of the ith molecule, expressed as 






(1) 


In this expression, oti is the permanent molecular polarizability of the ith molecule in isola¬ 
tion in the laboratory frame, and Tij is the dipole-dipole interaction tensor 


'■o 


1 ^ rij ® Tij 

o c 

T-- T-- 


( 2 ) 


Here, Vij is the vector from the center of mass of molecule i to the center of mass of molecule 
j, and Tij = \rij\. Also, 1 and ® are the unit matrix and the tensor product, respectively. In 
order to properly take into account the effect of the long-range interaction on the molecular 
polarizability, we employ the Ewald sum for Tij. This effect has been ignored in previous 
MD simulations of 2D Raman and 2D THz-Raman spectroscopy systems. However, the 
contribution of this effect in the 2D THz-Raman case is minor in comparison with that in 
the 2D Raman case. 

The total dipole moment is evaluated as ii(t) = where and 

are the permanent and induced molecular dipole moments of molecule i, respectively. 
The induced molecular dipole moment is expressed in terms of the interaction tensor as 


= Ri 


E. 


perm 






(3) 


where is the electrostatic held at molecule i created by all the other molecules in the 

system. This is evaluated as where Vu. is the vector between 

the center of mass of molecule i and that of atom I in molecule j . 


C. One- and two-dimensional signals 

It should be noted that the majority of MD simulations of 2D IR spectroscopy systems 
performed to this time have been carried out to obtain the parameter values for stochastic 
models. Full MD simulations have been carried out mostly in the cases of low frequency 
vibrational modes.®®' This is because the primary inter-molecular modes, which are the 
objects of study in 2D IR spectroscopy, can be separated from the other modes rather easily, 
as in the case of the OH stretching mode in liquid water. Contrastingly in the 2D THz- 
Raman case, it is not easy to hnd primary modes, because the objects of study in this 
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cases are inter-molecular vibrations that depend on the complicated nature of molecular 
ensembles, whose conhgurations change in time. Thus, we have to evaluate 2D signals 
directly from the MD simulations. Because quantum mechanical effects are minor for low- 
frequency inter-molecular modes, due to their small thermal activation energies, unlike the 
case of intra-molecular motion,®' and because 2D THz-Raman spectroscopy employs the 
three-body correlation function with two time variables, instead of the four-body correlation 
function with three time variables employed in 2D IR spectroscopy, the full MD simulation 
approach is practical. For this reason with our approach, we were able to carry out full MD 
simulations to directly evaluate 2D signals. 

Although our MD and model calculations are fully classical, we start from the quantum 
expressions for the response functions, because their classical expressions in the MD and 
model calculations are most easily derived by taking the classical limit of the quantum ex¬ 
pressions. The optical observables in ID and 2D spectroscopies are represented respectively 
expressed by two- and three-body response functions of the formJ^'^ 

R{t) = l{[A{t),m]) ( 4 ) 


and 

([[i(i2 + «i),S(ii)],C(0)]), (5) 

where A, B, and C can be the total dipole moment, fi, or the total polarizability of the 
molecules, 11. For low-frequency inter-molecular vibrations, we can take the classical limit, 
h —?• 0. The commutator and operators are then replaced by the Poisson bracket and 
c-number observables as 


-^[A,B]^{A,S}pb 


dA dB dA dB 
dq dp dp dq ' 


( 6 ) 


Using {e A(t)}pB = /9e for a molecular Hamiltonian HQ{p,q), we ob¬ 

tain the expression for thej linear response function, for example, for ID IR as 


^Ir(^) = /5(Req(t)Aeq(0)) 


(7) 


and 



(cu) oc culm 





( 8 ) 
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It should be noted that the quantum correction factor, tanh(/dha;/2), is usually included in 
Eq. ([^ to allow comparison with experimentally obtained IR signals, but here we do not 
include it. Instead, as the classical limit, we only multiply by We can evaluate the 
above quantities easily by calculating ^eq(g(t)) from samples of molecular trajectories q{t) 
that are obtained from the equilibrium MD simulation. 

For the 2D case, the response function in the classical limit is expressed ad^^® 

Rihiti) = ({{A(f2), .B(0 )}pb, C(—fi)}pB) (9) 

As in the ID case, we can calculate the above response function using iJ,{t) and n(f) eval¬ 
uated from the molecular trajectories p{t) and q{t) obtained from the equilibrium MD 
simulations.^I®^ The convergence of the signal is, however, very slow due to the effect of 
the stability matrix element in the double Poisson brackets. However, there is different ap¬ 
proach, the non-equilibrium hnite held approach, that does not have this convergence prob¬ 
lem. In this approach, the double Poisson brackets are evaluated as —{A(t'),B(t)}pB = 
(A+_B(t)(^0 “ where A±B{t)(t') is the observable corresponding to A(t') cal¬ 

culated from the trajectories subjected to the weak perturbations ±(—F5(r — t)S(r)), with 
the electric held ±F acting on .B(r).f23I25| approach is computationally intensive, 

and for this reason, here we employed a hybrid approach, which utilizes both the equilibrium 
and non-equilibrium approaches in order to reduce the computational cost further.^ 

In our hybrid approach, we evaluate C{—ti) = dC(t)/dt\t=-ti with equilibrium MD 
simulations, while A±B{o){t 2 ) with non-equilibrium MD simulations. As a result, the hybrid 
expressions for the 2D Raman-THz-THz, THz-Raman-THz, and THz-THz-Raman signals 
become 

-^RTt(^25H) = ^((r+/x( 0)(^2) ~ R-/.t(0)(^2)) rieq(—fi)), (10) 

3 2/9 

Rplj:{t2,tl) = -^r^((R+n(0)(^2) - R-n{0)(^2)) Aeq(-tl)), (11) 

and 

RTTR(^2,tl) = ^((n+^(o)(t2) - n_^(o)(t2)) Aeq(-H)), (12) 

where Ej and f3 = l/k-QT are the external electric held of the jth pulse and the inverse 
temperature divided by the Boltzmann constant. 







III. MODEL CALCULATION 

A. Brownian oscillator model with nonlinear interaction 


In order to analyze the 2D signals of inter-vibrational modes, we consider a model that 
consists of a primary oscillator mode nonlinearly coupled to the other modes, which are 
regarded as a bath system. This bath system is represented by an ensemble of harmonic 
oscillators. The primary mode may change in time or be inhomogeneously distributed. 
We can describe both situations within a unihed framework by adjusting the bath param¬ 
eter variables. The model is constructed by extending a Brownian (or Caldeira-Leggett) 
HamiltoniaiP^^ to include a nonlinear system-bath interaction. We write 

^ + Hi, (13) 


where 


is the Hamiltonian for the system with mass m, momentum p and potential U{q), 

^2 




2m 


1 




2mju‘j 


(14) 


(15) 


j ' ' j 

is the bath Hamiltonian with the momentum, coordinate, mass, and frequency of the jth 
bath oscillator given by pj, Xj, mj and Uj, respectively, and 


Hi =-V{q)'^ajXj, (16) 

j 

is the system-bath interaction, which consists of linear-linear (LL) and square-linear (SL) 
system-bath interactions, V{q) = 14l9 +Vsl?^/ 2, with coupling strengths Wl, Usl, and 
This model has been used to derive predictions for 2D RamarP^HSZ! ^nd 2D IR signals 
The last term of the bath Hamiltonian is the counter-term, which maintains the translational 
symmetry of the system in the case U{q) = 0. 

The sum of the bath coordinates X = ^2^ ctjXj acts as a collective coordinate that 
modulates the system.!^ As illustrated in Ref. EZl while the LL interaction shifts the po¬ 
tential, the SL interaction changes its curvature. Although in the anharmonic potential 
case, the LL interaction also changes the curvature of the potential, we can ignore this ef¬ 
fect if the anharmonicity is weak.^^^I^Sl Next we introduce the spectral distribution function. 
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(a) Fast modulation limit 
(Homogeneous case) 


(b) Slow modulation limit 
(Inhomogeneous case) 






FIG. 1. Schematic illustration of the relation between line broadening and modulation of the 
potential system perturbed by the SL system-bath interaction. (a)The fast modulation limit cor¬ 
responds to the homogeneous broadening case, whereas (b) the slow modulation limit corresponds 
to the inhomogeneous broadening case. 


J{uj) = — Uj)/2mjUJj, which characterizes the bath and system-bath coupling. 

We assume that J{uj) has an Ohmic form with a Lorentzian cutoff®®' 


J{u) 


hm( 

2 tt 7^ -f ’ 


(17) 


where C is the system-bath coupling strength, and 7 represents the width of the spectral 
distributuion. 

Writing the classical collective coordinate corresponds as W, we have the correlation 
function, (df(t)X(O)) oc This indicates that the bath oscillators interact with the 

system in the form of Gaussian Markovian noise with correlation time r = I/ 7 .® Because 
the SL interaction affects the frequency of the potential, the fast modulation limit (r —>■ 0) 
of the SL interaction corresponds to the case of a homogeneous distribution, as depicted in 
Fig. 0 a). By contrast, the slow modulation limit (r —)■ 00 ) corresponds to the case of an 
inhomogeneous distribution, as depicted in Fig. [^b). Note that SL or LL-I-SL BO model 
has the “mode mixing in polarization”® because the modes included in collective coordinate 
q are frequency distributed by SL interaction and are mixed by nonlinear polarizability 1120 '^. 
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B. Classical hierarchal Fokker-Planck equations 


Because we wish to explore the effects of anharmonicity, nonlinear polarizability, vibra¬ 
tional dephasing, and homogeneous and inhomogeneous broadening within a unified frame¬ 
work, we must employ a kinetic equation that can treat thermal fluctuations as well as 
dissipation in a non-perturbative, non-Markovian manner. The reduced hierarchal equa¬ 
tions of motion (HEOM) satisfy all of the requirements mentioned above and is ideal for the 
present study.®® While we must use the quantum form of the equations to calculate the 
signals for high frequency intra-molecular modes,® we can employ the classical form for low 
frequency inter-molecular modes. 


For the LL-I-SL BO Hamiltonian, given in Eqs. (13)-(16), the HEOM for the classical 
distribution function are expressed a^®® 


dt 


= - (^£ci + ny) (p, q; t) - ny©(p, g; t) - (p, q; t) (18) 


for 0 < n < iV and 


dW(^\p,q;t) ^ _ NjQW^^-^\p,q-t). (19) 

In the HEOM approach, only the first element, W^^\p,q;t), has physical meaning, while 
the other elements, W^'^\p,q;t) (1 < n < N), are introduced in the numerical calculations 
in order to treat the non-perturbative, non-Markovian system-bath interaction. We choose 
N to satisfy N 3> ujc/l, where Uc is the characteristic frequency of the system. 

The classical Liouvillian of the system, Cci, is defined by 


m oq op 


( 20 ) 


where the dash is defined as A'{q) = dA{q)/dq for an arbitrary function 24(g). The operators 
<h and 0 describe the energy exchange between the system and the heat bath for the inverse 
correlation time y. They are defined as 


and 




0 ^ -CV(,) (p + 


( 21 ) 


11 


( 22 ) 






with /3 = l/k-^T and 


V\q) = liL + Vsi^q. 


(23) 


The thermal equilibrium distribution, q), is expressed in terms of the HEOM elements 


evaluated from the steady-state solution of the HEOM. Note that Eqs. (18)-(19) reduce to 
the Kramers equation in the limit iV —)■ 0 with Vsl = 0.^ 

Hereafter, we employ the dimensionless coordinate and momentum defined by g = 
uo^/ml]/2 X q and p = x p, where ooq = \JU''{q)/m represents the fundamen¬ 

tal frequency. The potential is then assumed to be 


U{q) = + ^9’, 


(24) 


where is the cubic anharmonicity of the potential. The other variables, Ell , Vsl, p{q), 
and n(g), are also normalized accordingly. 

C. One- and two-dimensional signals 


To apply the HEOM formalism, we express the response functions in terms of the time- 
propagation operator. Then, Eqs. (|^ and ([^ can be rewritten as 

(25) 


R(t) = -Tt 


and 


R{t2,ti) = ( ^ ) Tr 




(26) 


where we have employed the hyperoperator ^ dehned as A^B = Q(t) is the Green’s 

function of the system Hamiltonian without a laser interaction, and peq is the equilibrium 
state. The above equations represent the time evolution of the system under laser excitation. 


For example, Eq. (26) can be interpreted as follows. The system is initially in the equilibrium 


state and is then modihed as a result of the hrst laser pulse via the dipole interaction 
by C. It then propagates for time ti under Q{ti). The system is next excited through the 
second laser pulse by B and propagates for time t 2 under ^(^ 2 )- Finally, the expectation 
value of the polarizability at H -|- ^2 is generated through the laser pulses by AW^ 
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The classical expressions for the response functions can be obtained from the above with 
the use of the Wigner transformation.^ In this case, an arbitrary operator A^{q) is replaced 
by A'{q){d/dp). For ID IR and 2D THz-Raman spectroscopies, we have 


^iR (t) = I I dgg 


and 


R 


r,Tt(^ 2, ^l) — h-ilfi 


d 




^TRT(^2,il) = piAi 


fipj 

f dgg 

[dpj 

f dgg 


d 




d 


Gih) {l + n2q)^{Gih)-W^\p,q) 


dp 

d 


and 


^ttr(^2,^i) = pfAi / dp / dg ( g + 


d 


dp 


dp 


d 


Git2)^[Giti)^W^\p,q) 


dp 


(27) 

(28) 

(29) 

(30) 


respectively. Here, II 2 = n 2 /ni is the relative intensity of the nonlinear polarizabil¬ 
ity. The Green’s function G{t) is now expressed in terms of the classical Liouvillian as 


G{t) = exp 


-Celt 


, and q) is the equilibrium distribution. In order to apply 

the HEOM formalism, we express the time-dependent Wigner function W{p,q-,t), such as 
^{p^qAi) = G{ti)dW‘^‘^{p,q)/dp and W{p,q;ti + t 2 ) = Git 2 )dW{p, q;ti)/dp, in terms of 
the HEOM member g; ti) and g;ti -|- ^ 2 ), and the determine its time evolu¬ 
tion through Eqs.(18) and In the HEOM formalism, the equilibrium distribution, 

g), is also expressed using the HEOM elements evaluated from the steady-state so¬ 


lution of Eqs. (18) and (19). In the strong coupling case, we employ an eigenfunctional 
representation of the momentum space for numerical convenience, as discussed in Appendix 
[A| while, in other cases, we solve Eqs. (18) and (19). Note that, as shown in Ref. [3ll 
the hierarchal Fokker-Planck approach is equivalent to the generalized Langevin approach. 
However, because the Fokker-Planck approach does not require sampling of system trajecto¬ 
ries, unlike the Langevin approach, it is numerically advantageous, especially for calculating 
nonlinear response functions, for which the trajectories are unstable. 


IV. RESULTS AND DISCUSSION 

Because the 2D profiles of the signal must be constructed from complex motion in a 
complicated manner, the analysis of the signal profile is not straightforward. Nevertheless, 
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properly accounting for the components of signals to these prohles allows us to perform a 
detailed analysis of the inter-molecular vibrational motion on the basis of both experiential 
and theoretical results. Model-based studies of the 2D prohles of signals are helpful to 
identify the underlying physical mechanisms, because it is necessary to capture the essential 
features of the inter-molecular motion in order to reproduce the complex 2D prohle from 
a simple model. Analyses of this kind have employed LL -|- SL BO models for the 2D 
Ip[22H29|49] 2D Raman cases.®® However, their applicability has not been fully explored 
because of the limited availability of experimental and numerical data. 2D THz-Raman 
measurements, which are applicable not only to the study of liquids, but also to the problem 
of distinguishing optical processes through use of RTT, TRT, and TTR measurements, 
provide the opportunity to explore the possibilities of model-based analysis of 2D prohles. 
Here, we examine the LL-I-SL BO model and use it to reproduce all three THz-Raman signal 
prohles obtained from full MD simulations. This is done by choosing the parameter values 
of the model so as to realize the best agreement between the signal prohles provided by the 
model and the MD simulation. Before htting the model to the MD signals, we demonstrate 
a general aspect of 2D THz-Raman signals using the SL BO model and the optical Liouville 
paths in a simple case. This serves as a guid to subsequent analysis. 


A. General aspects of 2D THz-Raman signals 


Here, we elucidate several aspects of the signal components in 2D THz-Raman spec¬ 
troscopy using the SL BO model and optical Liouville paths. As explained in Sec. we can 
express the dipole and polarizability in terms of the collective coordinate q as i^{q) = fiiQ 
and n(g) = U,q + n2g^/2, respectively. If the inter-molecular modes are both Raman and 
THz active, the three THz-Raman signals are expressed as 


and 



p(3) 

n-TTR 


{h, ti) 


njUi 

fijUi 

/iiHi 


Rah(0, H) + 
-Rah( 0 , H) + 


Rah( 0 , H) + 


-^Rrtt(0,H) , 
-^-Rtrt(0,H) , 


fh 

2 


Rttr(0,H) , 


(31) 

(32) 


(33) 
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with 112 = n 2 /ni, where the a n h arm on in component is expressed as 

RAK{t2, = {[[q{ti2),q{ti)], g(0)]). (34) 

Similarly the nonlinear polarizability components are given by 

RKrT{t 2 ,ti) = -^ {[[q{t 12 ), q{ti)],q'^{0)]) , (35) 

-Rtrt(^ 2 ,^i) = ([[ 9 (^ 12 ), g^(^i)],g(0)]) , (36) 

and 

R-TTRih, ti) = —^ ([[9^(^i2), 9 (^ 1 )]) 9 ( 0 )]) • (37) 

In the harmonic LL BO case, the above THz-Raman signals can be calculated analytically, 
and we have Rah = 0, Rrtt = 0, 

Rtrt(ORi) = —^CRRCih), (38) 

and 

Rttr(05^i) = —^C{ti + t2)C{t2), (39) 

where C{t) = {[q{t), q]) is the hrst-order response function of the harmonic BO system. For 
an isolated oscillator with frequency u, we have C{t) = hsm{ut)/2mu. As explained in 
Sec. 1^ the term Rah vanishes in the harmonic case because there is an odd number of 
Gaussian integrals involved in the response function.l^^^^ Moreover, the term Rrtt vanishes 
because of the cancelation of possible optical Liouville paths, as explained in Appendix 
While Rah becomes large for large anharmonicity the contribution of Rrtt remains 

small due to this cancelation. Thus, in the anharmonic case, we can estimate Rah from the 
RTT measurement. Then, by subtracting Rah from the TRT and TTR signals, Rtrt 
Rttri evaluate Rtrt and Rttr separately. Because each contribution arises from 

the corresponding optical process, we can elucidate the key features of the inter-molecular 
motion from them. In the 2D Raman case, contrastingly, because all of the contributions 
appear together in a single observable, we cannot carry out such analysis. 

To more clearly elucidate the characteristic of the 2D THz-Raman signals, in the hgures, 
we display the 2D prohles of the Rtrt, Rttr, and Rah components separately in the slow 
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TABLE I. Relative intensities of the signal components for the SL BO model in the harmonic 
(53 = 0) and anharmonic (g^ = 0.3) cases. The intensities are estimated from the maximum peak 
value of the signal normalized with respect to the values of Rtrt- Although the SL interaction 
gives rise to the contribution of Rrtt to the RTT signal, its intensity is weaker than that of Rah 
if we include the prefactor II 2 . This situation does not change even if we increase the strength of 
the anharmonicity. 


Potential 

Modulation 

Rah Rrtt“^ 

Rtrt Rttr 

harmonic 

slow 

— 

0.18 

1 

0.99 

harmonic 

fast 

— 

0.11 

0.85 

0.81 

anharmonic 

slow 

0.12 

0.18 

1.01 

1.00 

anharmonic 

fast 

0.14 

0.13 

0.85 

0.81 


The 2D profile of the RTT component is presented in Appendix [b] 

modulation case = 1.0 Uq and 7 = 0.5 Wq) and the fast modulation case = 0.49 ojq and 
7 = cx)), as obtained from the SL BO model (Vll = 0, Vsl = 1) with ujq = 600 cm“^ and 
T = 300 K. Note that, because the effective system-bath coupling strength becomes weaker 
as the modulation becomes faster, we change both 7 and C, to elucidate a pure non-Markovian 
effect.^ Although the contribution is minor, we plot the Rrtt component in Appendix 
The relative intensities of each component evaluated in the harmonic and anharmonic SL 
BO cases are presented in Table. We then analyze each prohle using the optical Liouville 
paths (the double-sided Feynman diagrams). 


1. The TRT component, i?TRT(05^i 


In Fig. 1^ we present the Rtrt component for (a) the slow and (b) the fast modulation 
cases calculated using the harmonic SL BO model {g^ = 0). It is seen that the peak prohles 
are symmetric along the ti = t 2 line, as can be deduced from the analytical expression for 


the LL case, Eq. (38). The peaks in the slow modulation case stretch in the — ^2 = 2 n 7 r/a; 


direction, whereas those in the fast modulation case stretch in the - 1-^2 = 2 n 7 r/a; direction, 
where u is the fundamental frequency and n an integer. In the slow modulation case, there 


16 










FIG. 2. The 2D profiles of the .Rtrt component calculated using the harmonic SL BO model 
(<73 = 0) for the (a) slow and (b) fast modulation cases. Contours in red and blue represent 
positive and negative values, respectively. 
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FIG. 3. The Liouville paths involved in .Rtrt- The red circles and the blue double circles rep¬ 
resent single and double quantum transitions, respectively. The three other paths, which are the 
Hermitian conjugates of the above-mentioned paths, are not presented. These conjugate paths can 
be obtained by exchanging the left and right arrows. 


are elongated peaks called ’’echo peaks” along the ti = t 2 direction. 

Althongh our simulation results are fully classical, these prohles can be interpreted easily 
using the quantum Liouville paths for the TRT process depicted in Fig. There, an energy 
eigenstate of the harmonic potential is denoted by |n), and we have depicted cases starting 
from the vibrational ground state as examples. The dipole operator pig, represented by the 
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red circles in the diagram, converts the state \n) into \n + 1) and \n — 1) through single 
quantum (SQ) excitations, while the nonlinear polarizability operator 112^^, represented by 
the blue double circles, converts the state \n) into |n+2) and \n — 2) through double quantum 
(DQ) excitations or maintains the same state |n) through the zero quantum (ZQ) excitation. 
Because the hnal state, appearing after the last laser interaction, must be a population state 
|n){n|, due to the trace operation involved in the response function, the possible processes 
are limited to the cases depicted in Figs, [^i)-(iii) and their conjugate diagrams. While the 
double circle in the diagram (i) involves the transition {a^a + aa^)|n + 1) = {2n + 3)|n + 1), 
that in the diagram (ii) involves the transition (a^a + aa^)|n) = {2n + l)|n). Thus, although 
these paths have same phases with opposite signs, they do not cancel. These six components 


constitute the signal expressed in Eq. (38) in the isolated oscillator case. 

The phases for the paths (i) and (ii) in Fig. |^i), expressed as exp[—za;(fi + ^ 2 )], are the 
same because the net transition for these paths is a ZQ transition, while that for path (hi), 
expressed as exp[—ia;(ti — ^ 2 )], is different because the net transition for this path is a DQ 
transition. In the slow modulation case, the signal in Fig. (iii) can rephase and become 
strong along the ti = ^2 direction, but in the fast modulation case, coherence is lost due to 
rapid changes in the fundamental frequency, as illustrated in Fig. and the signal decays 
quickly. Thus, we observe the echo signal in the TRT component in the slow modulation 
case, whereas we observe the chain of peaks along the ti + t 2 = 2n7i/u direction in the fast 
modulation case, due to the processes depicted in Figs. (i) and (ii). 


2. The TTR component, i?TTR(^ 2 ,^i) 


We now discuss the Rttr component, presented in Fig. The prominent features of 
this signal are the appearance of elongated peaks parallel to the ^2 axis and the appearance 
of peaks along the ti + 2^2 = 2n'n ju direction for small values of ^ 2 - The optical Liouville 
paths for Rttr are presented in Fig. Three other paths can be obtained by exchanging 


the left and right arrows. These six components constitute the signal expressed in Eq. (39) 
in the isolated oscillator case. 

We observe the population states \n){n\ during the ^2 period of Figs. |^i) and|^(ii) and 
the coherent states |n + 2)(n| created by the two SQ excitations during the t 2 period of Figs, 
[^iii). While all of the diagrams in Fig. [^exhibit the oscillation exp[—icuti] during the ti 
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FIG. 4. The 2D profiles of the .Rttr component calculated using the harmonic SL BO model 
(<73 = 0) for the (a) slow and (b) fast modulation cases. Contours in red and blue represent 
positive and negative values, respectively. 



FIG. 5. The Liouville paths involved in .Rttr- The red circles and the blue double circles rep¬ 
resent single and double quantum transitions, respectively. The three other paths, which are the 
Hermitian conjugates of the above paths, are not presented. 


period, only that in Fig. |^iii) exhibits the oscillation exp[—2zci;t2] during the t 2 period. Due 
to the SL interaction, the high-frequency oscillation of the coherent state appearing in Fig. 
[^iii) decays quickly, while the population state appearing in Fig. [^ii) remains for a long 
time. Thus, we observe the peaks along the ti + 2^2 = ^rni/uj direction for a short t 2 period, 
whereas elongated peaks appear along the t 2 direction. 
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FIG. 6 . The 2D profiles of the .Rah component calculated using the anharmonic SL BO model 
(<73 = 0.3) for the (a) slow and (b) fast modulation cases. Contours in red and blue represent 
positive and negative values, respectively. 


3. The AH component, Rah(^ 2 ,^i) 

The contribution of the component Rah arises only in the anharmonic case. In Fig. 
we display the signal for = 0.3 in the (a) slow and (b) fast modulation cases, respectively. 
Although the definition of the response function is different, the optical Liouville paths for 
Rah are similar to those for Rtrt and Rttr- To explain the reason for this, we consider 
energy eigenstates in the case of an anharmonic potential In') with eigenenergy Un'- The 
fundamental frequency between | 0 ') and | 1 ') is denoted by u' = uy — cjo', whereas that 
between |1') and |2') is denoted hy uj' — 5 = U2' — ojy. In addition to such frequency shifts, 
the anharmonicity changes the roles of the dipole and polarizability operators. While the 
Rah component involves only piq and Ilig, they can induce ZQ and DQ transitions, in 
addition to the SQ transition, because {n'\q\n'), {n'\q\n' + 2 '), and {n' + 2'|g|n') are all 
nonzero in the anharmonic case. Because both THz and Raman laser pulses can induce ZQ 
and DQ transitions, the diagrams for Rah include all of the diagrams in Figs. [^and[^ while 
the resonant frequency for the ti period becomes ca', and that for the t2 period is 0 , u', 
oj' — 6, or 2a;' — 6. The rephasing processes depicted in Fig. (iii) occur only rarely, even 
in the slow modulation case, due to the anharmonicity, and this contribution can therefore 
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TABLE II. The fitted parameter values of the LL+SL BO model for formamide, water and 
methanol liquid. Here, have fixed = 0.1, gi = 1, Hi = 1, and VsL = 1- 


Molecule 

iOo 

(cm 

n2 

C/<^o Fll/Fsl 

Formamide 


60 

0.04 

20 

oo 

0 

Water 


450 

0.07 

6.0 

0.7 

0.01 

Methanol 


540 

0.07 

2.1 

0.4 

0.01 


be ignored. For this reason, the 2D prohles of the Rxr component are similar to those of 
-Rttr and .Rttr, bnt without the echo peaks described by Fig. (hi). As in the TTR case, 
the contribution from the diagram depicted in Fig. (hi) decays quickly. Thus, we observe 
elongated peaks in the ^2 direction that arise from population decay during the t 2 period. 


B. The MD and fitted results 

In order to concretely study hydrogen-bonding dynamics in 2D THz-Raman spectroscopy, 
we selected formamide, water, and methanol. While water and methanol exhibit high- 
frequency librational motion arising from hydrogen bonding, formamide exhibits only low- 
frequency inter-molecular motion. For each liquid, we explore the parameter values of the 
LL-I-SL BO model to £t the RTT, TRT and TTR signals to the MD simulation results. We 
hxed the anharmonicity to = 0.1, while we varied the nonlinear polarizability 112. This 
was done because the intensity of the Rah component is proportional to and only 

the ratio of and n 2 is important to elucidate the roles of anharmonicity and non-linear 
polarizability. Because the coupling strengths of the LL and SL interactions are determined 
by C^L and C^sl, respectively, we also fixed Vsl = 1 and varied C and VLu/hsL in fhe 
£t. The best hts of these parameters are presented in Table |TTj Note that because the SL 
interaction increases the effective system frequency,®® the fitted uq do not correspond to 
the resonant frequency estimated from the ID IR results obtained from the MD simulation. 
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Formamide 


(a) MD simulation 


(i) RTT (ii) TRT (iii) TTR 



(b) BO model 





(i) RTT 


(ii) TRT (iii) TTR 



h (ts) 


(c) Linear absorption (IR) 



FIG. 7. The 2D signals for formamide obtained from (a) the MD simulation and (b) the LL+SL BO 
model for the (i) RTT, (ii) TRT, and (iii) TTR measurements, respectively. The fitted parameter 
values of the LL+SL BO model are presented in Table |IIj Using these values, we also calculated 
the ID IR signal, which is presented with the MD results in (c). 
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1. Formamide 


In Fig. we display (a) the MD results and (b) the htted results from the LL+SL BO 
model for the (i) RTT, (ii) TRT, and (hi) TTR cases of liquid formamide. It is seen that 
the LL+SL BO results capture the essential features of the 2D THz-Raman spectra. The 
characteristic feature of the MD signals for formamide is the elongation of the peaks along 
the t 2 axis observed in Figs, [^a-i) and (a-ii). By adapting the analysis applied to Rah and 
Rtrt, we conclude that this elongation arises from the slow population decay during the 0 
period described by the diagrams in Figs. (ii), 1^ (i), and[^ (ii). The echo peaks do 

not appear in Figs, [^a-ii) and (b-ii), because the modulation in this case is so fast that 
dephasing cannot take place. In the TTR case depicted in Figs, [^a-iii) and (b-iii), the 
peaks along the ti + 2^2 = 2n7r/u direction appear for small values of ti and ^2 due to the 
contribution of the diagram presented in Fig. (iii). The similarity between the MD and 
model results indicates that the collective modes are subject to fast frequency modulation, 
as illustrated in Fig. [^a). 


2. Water and Methanol 

In Figs. [^and|^ we display (a) the MD results and (b) the LL+SL BO results for the 
(i) RTT, (ii) TRT, and (iii) TTR cases of liquid water and methanol, respectively. We 
reproduced all of the 2D prohles obtained from the MD simulation for these liquids using 
the LL+SL BO model with the parameter values listed in Table In these liquids, the 
vibrational modes near oj = 600 cm“^ are the librational modes.^^^® 

The characteristic features of the MD signals for these two liquids are the elongation of 
the peaks along the ^2 axis in the RTT and TRT signals and the echo peaks along the ti = t 2 
line in the TRT signals. The existence of these peaks indicates that the collective modes are 
inhomogeneously distributed, as depicted in Fig. (b), due to slow frequency modulation. 
Because water and methanol exhibit high-frequency librational motion caused by hydrogen 
bonds, we observe an oscillatory feature in the 2D signal. This contrasts with the formamide 
signal, which decays quickly. 
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Water 


(a) MD simulation 


(i) RTT (ii) TRT (iii) TTR 



(b) BO model 

(i) RTT (ii) TRT (iii) TTR 



(c) Linear absorption (IR) 



0 200 400 600 800 1000 1200 


CO (cm ^) 


FIG. 8. The 2D signal for water obtained from (a) the MD simulation and (b) the LL+SL BO 
model for the (i) RTT, (ii) TRT, and (iii) TTR measurements, respectively. The fitted parameter 
values of the LL+SL BO model are presented in Table |IIj Using these values, we also calculated 
the ID IR signal, which is presented with the MD results in (c). 
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Methanol 


(a) MD simulation 

(i) RTT (ii) TRT (iii) TTR 



0 100 200 0 100 200 0 100 200 


h (fs) 


(b) BO model 

(i) RTT 


(ii) TRT 


(iii) TTR 



1-r 


0 




0 100 200 0 


100 200 0 
h (fs) 


100 200 


(c) Linear absorption (IR) 



FIG. 9. The 2D signal for methanol obtained from (a) the MD simulation and (b) the LL+SL BO 
model for the (i) RTT, (ii) TRT, and (iii) TTR measurements, respectively. The fitted parameter 
values of the LL+SL BO model are presented in Table |IIj Using these values, we also calculated 
the ID IR signal, which is presented with the MD results in (c). 
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V. CONCLUSION 


Using a Brownian oscillator (BO) model with a linear-linear (LL) and square-linear (SL) 
system-bath interaction, we analyzed the 2D RTT, TRT, and TTR signals for formamide, 
water, and methanol liquids obtained from full MD simulations. The classical hierarchal 
equations of motion approach was used to calculate the 2D signals with the LL-I-SL BO 
model under non-perturbative and non-Markovian conditions. By htting the anharmonicity 
of the potential, the LL and SL coupling strength, and the inverse noise correlation time 
of the LL-I-SL BO model to the results of the MD simulations, we reproduced each of 
the simulated signal prohles by capturing their characteristic features. We found that the 
prohle of formamide liquid can be accounted for by homogeneously distributed oscillators, 
whereas the prohles of water and methanol liquids can be accounted for by inhomogeneously 
distributed oscillators. Due to hydrogen bond interactions, water and methanol exhibit 
oscillating echo peaks. We were able to describe the two cases as the cases of fast and slow 
modulation of the anharmonic LL-)-SL BO model. The key feature of the present model 
that allows it to account for the simulated signal is the existence of the SL interaction. We 
were able to describe the complex inter-molecular motion with this simple model because 
it is capable of describing the collective modes under conditions ranging from homogeneous 
to inhomogeneous in a unihed manner through variation of the noise correlation time. The 
success of the present study indicates that the LL-I-SL BO model captures the essence of 
the inter-molecular motion. 

Finally, we briefly discuss some extensions of the present study. First, note that there 
does exist some discrepancy between the MD simulation results and the LL-I-SL BO model 
results. It should not be too difficult to decrease this discrepancy. For example, in the case of 
water, we should be able to improve the description of the signal prohles by introducing the 
second mode, which is 220 cm“^, in the BO model. The signihcance of the htted frequencies, 
(Uq, should be clarihed through the normal mode analysis of MD simulations. Secondly, as 
shown in previous studies, the LL-I-SL BO model can be applied not only to 2D THz-Raman 
spectroscopy, but also to 2D IR spectroscopy for both inter-moleculai^ and intramolecular 
vibrations. extension to bridge between the inter- and intra-molecular modes using 
an extension of the present model should also be possible. We leave such extensions to future 
studies, in accordance with progress in experiments and simulations. 
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Appendix A: The Hierarchal Fokker-Planck Equations in the Hermite 
Representation 


In the case of a strong system-bath coupling strength Eqs. (18) and (19) converge very 
slowly as difference equations with a discrete mesh in the phase space. By expanding in 
terms of Hermite functions in the momentum direction,^ the convergence can be improved. 
In this case, the HEOM become simultaneous equations for the coefficients of the expression. 
We expand the distribution function as follows: 


W^'^\q,p]t) =^IJoe 


where ‘ipk{p) is the kth Hermite function, 

= 


k=0 


p 


=Hk ( - ) exp 


p 


2a2 


(Al) 


(A2) 


with Hk{x) is the kth Hermite polynomial and a = \j2m/{3. 

The Liouvillian of the system (given in Eq. (20)) and the relaxation operators (given in 


Eqs. (21) and (22)) are expressed as 


L = + 6+D~, 

$ = V'(q)\ —b+, 

V m 


0 = 



where A = ^3p^A^+3Ul2^^-i3p^/4m-i3ui2^ 


_ 1 n3 fm d 

2 V ra^ ^ V (3 dp'' 


(A3) 

(A4) 


(AS) 


(A6) 
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and 


_ 1 hS dU 1 jm d 
2\mdq m\J (3 dq 


(A7) 


Then, the eqnations of motion for the coefficients c^i^\q ■ t) are rednced to 


dc 


(n) 


(n) 


dt 


Vk + lD^c^kii - VkD c^k-i - 




(A8) 


for 0 < n < 77 and 


dc 


{N) 


^ ^Jk + 10+47 - 'fio 4-1 - ^7C, 


(«) 


+ iV7CV"(«)7^7F + T47‘> - (,V'(qf-kA^'- 


(A9) 


To carry ont the nnmerical calcnlation, we chose fcmax so as to satisfy ~ 0 (fc > fcmax) and 


solved Eqs. (A8)-(A9) as /cmax simnltaneons eqnations. 


Appendix B: The RTT component, i?RTT(^25^i) 


In Fig. 10, we depict the Rrtt component for the (a) slow and (b) fast modnlation cases, 
calcnlated nsing the harmonic SL BO model {g^ = 0). Here, we chose the same parameter 
valnes as in IV-A. The estimated signal intensity for Rrtt is presented in Table 

First we shonld note that the signal intensity of this component is very weak. To illnstrate 


this point, we present the optical Lionville paths for Rrtt in Fig. 11 In the harmonic case. 


the diagrams in Figs. 11'i) and 11 (ii), and those in Figs. 11 qn) and 11 (iv) 

cancel, respectively, and the signal from Rrtt vanishes. However, in the anharmonic case 


and/or the SL BO case, the diagrams in Figs. [lT[iii) and 11 (iv) may snrvive, althongh those 
in Figs. BHi) and (ii) always cancel. This is becanse the transition freqnencies between |0) 


and |1) and between |1) and |2) involved in the ^2 period of Figs. [lT[iii) and 11 (iv) are 


different in these cases. Tims, we observe the signal displayed in Fig. [^with nodes along 
the 2ti + t 2 = rniju directions. Its intensity is weak, however, becanse of the cancelation. 
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FIG. 10. The 2D profiles of the .Rrtt component calculated using the harmonic SL BO model 
(<73 = 0) for (a) tire slow modulation case ((" = 1-0 ujq and 7 = 0.5 wq) and (b) the fast modulation 
case (C = 0.49 loq and 7 = 00 ). Contours in red and blue represent positive and negative values, 
respectively. 



FIG. 11. The Liouville paths involved in .Rrtt- The red circles and the blue double circles represent 
single and double quantum transitions, respectively. The four other paths, which are the Hermitian 
conjugates of the above paths, are not presented. In the harmonic case, the propagators from (i) 
and (ii), and those from (iii) and (iv) cancel, and the signal from this component vanishes. 
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TABLE III. Relative intensities of the signal components for the LL (Tll = Ij hsL = 0), SL 
(Ell = 0, VsL = 1), and LL+SL LL (Vll = 0.5, Vsl = 1) anharmonic BO model (<73 = 0.1). The 
intensities are estimated from the maximum peak values of the signals. The obtained intensities 
were normalized with respect to the intensity of Rtrt in the slow modulation case. 


Model Modulation 

Rah 

Rrtt 

Rtrt Rttr 

LL 

slow 

0.04 

< 0.01 

1 

1.00 

LL 

fast 

0.07 

0.01 

1.22 

1.22 

SL 

slow 

0.05 

0.23 

1.28 

1.26 

SL 

fast 

0.06 

0.14 

1.09 

1.04 

LL+SL 

slow 

0.37 

0.16 

1.15 

1.11 

LL+SL 

fast 

0.15 

0.10 

1.16 

1.13 


Appendix C: The 2D profiles of the signal components for different LL+SL 
couplings 


In this appendix, we present the profiles of the 2D THz-Raman signals in terms of the (i) 
AH, (ii) RTT, (iii) TRT and (iv) TTR components to clarify the roles of the LL and/or SL 
interactions. These results can be used to elucidate the key features of the inter-molecular 
interactions as revealed by 2D signals obtained from experiments or simulations. 


To carry out the numerical calculations, we set the frequency, anharmonicity, and tem¬ 
perature as cjq = 600 cm“^ , = 0.1, and T = 300 K. Then we calculated the signal 

components for the LL (1+l = 1, Rsl = 0), SL (1+l = 0, I+l = 1), and LL+SL (1+l = 0.5, 
Vsl = 1) models in the slow (7 = 0.5 wq) and fast (7 = cxd) modulation cases. Because the 
effective coupling strength depends on 7 , 1^21 we adjusted ( in the fast and slow modulation 
cases in order to have the same maximum intensity in the absorption spectrum. 


The 2D prohles of the (i) Rah, (h) Rrtt, (hi) Rtrt, and (iv) Rttr components for the 
LL, SL and LL+SL models are presented in Figs. 1^ ^ respectively. The estimated signal 
intensities for each component are listed in Table I 


111 


Because the anharmonicity is not strong, the calculated results appearing in Figs. 12 
andlT^iv) are similar to those predicted by Eqs. (38) and (39), respectively. This indicates 
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LL 


(a) slow modulation (b) fast modulation 


(iv) Rttr 
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FIG. 12. The 2D profiles of the (i) .Rah, (h) Rrtt, (hi) Rtrt, and (iv) Rttr components for the 
(a) slow modulation ((" = 1.5 ujq and 7 = 0.5 loq ) and fast modulation ((^ = 0.2 wq and 7 = 00) 
cases obtained from the anharmonic (g'3 = 0.1) LL BO model. Contours in red and blue represent 
positive and negative values, respectively. 
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(a) slow modulation (b) fast modulation 



(iii) Rtrt J 





FIG. 13. The 2D profiles of the (i) .Rah, (h) Rrtt, (hi) Rtrt, and (iv) Rttr components for the 
(a) slow (( = LOo and 7 = 0.5 wq ) and fast modulation ((" = 0.49 loq and 7 = 00) cases obtained 
from the anharmonic (53 = 0.1) SL BO model. Contours in red and blue represent positive and 
negative values, respectively. 
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(a) slow modulation (b) fast modulation 
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FIG. 14. The 2D profiles of the (i) Rar, (h) .Rrtt, (hi) .^trT) and (iv) .Rttr components for 
the (a) slow (C = 1-35 loq and 7 = 0.5 cjq) and fast modulation (C = 0.32 cao and 7 = 00) cases 
obtained from the anharmonic {g^ = 0.1) LL+SL BO model (Vll = 0.5 and Dsl = !)• Contours 
in red and blue represent positive and negative values, respectively. 
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that the node lines in -Rtrt correspond to ti = 2111 ^juj and = 2m7r/a;, whereas those 

in -Rttr correspond to ti + 2^2 = 2m^juj^ where n and m are any integers. Note that, 

due to the effect of the heat bath, the fundamental frequency uj is shifted from The 

components .Rah and .Rrtt in the LL case appear solely due to the anharmonicity, while 

those presented in Figs. and arise from both the SL interaction and the anharmonicity 
in the former case and just the SL interaction in the latter case. The profiles in the fast 
modulation limit of the SL mode presented in Figs, [^b) and 13|^b) are similar to those 
for the anharmonic LL model presented in Fig. because this limit corresponds to the 
case of homogeneously distributed oscillators, illustrated in Fig. [^a). The intensity of the 
Rrtt component in the LL case is significantly weaker than those of the other components. 


because of the cancelation explained in Appendix B. As in the case of Fig. we observe 
nodes along the 2t\ + ^2 = 2n'Kjuj direction, although they are not clear. 


Although we included anharmonicity in the SL case depicted in Fig. 13, the overall 
profiles of the signals are similar to those in the harmonic case, presented in Figs. 0 i and 


10 This indicates that the effects of the SL interactions are dominant in these profiles. This 


is also true in the LL+SL case considered in Fig. M, as indicated by the similarity of the 
profiles there with those in Figs. i0 [T0I and The intensity of the Rah component 
in the LL+SL case, however, is much larger than that in the other cases, because the 
LL+SL interaction induces transitions through the cubic interactions 

VLLVsLQ‘^it")Q(t'), derived from V{q(t"))V{q(t')). These interactions can induce a signal 
in the .Rah component, even in the harmonic case, involving only Gaussian integrals in 
the response function, for example, as lALTsLTr{q'(fi +t 2 )q{t")q{ti)q‘^{t')q{ 0 ) exp{—(]Hs)}^ 
Then, the profile of the .Rah component in the LL+SL case exhibits some similarity to 
that in the LL model, presented in Fig. [I^iii), because the dipole and linear polarizability 
operators associated with the system-bath interactions, for example, q(t")q(ti)q^(t') and 
q{ti + t 2 )q{t")q‘^{t'), induce the same DQ and ZQ transitions as nonlinear polarizability. 
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